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Abstract 

An important class of problems exhibits macroscopically smooth be- 
haviour in space and time, while only a microscopic evolution law is 
known. For such time-dependent multi-scale problems, the gap-tooth 
scheme has recently been proposed. The scheme approximates the evolu- 
tion of an unavailable (in closed form) macroscopic equation in a macro- 
scopic domain; it only uses appropriately initialized simulations of the 
available microscopic model in a number of small boxes. For some model 
problems, including numerical homogenization, the scheme is essentially 
equivalent to a finite difference scheme, provided we repeatedly impose 
appropriate algebraic constraints on the solution for each box. Here, we 
demonstrate that it is possible to obtain a convergent scheme without 
constraining the microscopic code, by introducing buffers that "shield" 
over relatively short times the dynamics inside each box from boundary 
effects. We explore and quantify the behavior of these schemes system- 
atically through the numerical computation of damping factors of the 
corresponding coarse time-stepper, for which no closed formula is avail- 
able. 

1 Introduction 

For an important class of multi-scale problems, a separation of scales exists 
between the (microscopic, detailed) level of description of the available model, 
and the (macroscopic, continuum) level at which one would like to observe 
the system. Consider, for example, a kinetic Monte Carlo model of bacterial 
chemotaxis iSCK03 | . A stochastic biased random walk model describes the 
probability of an individual bacterium to run or "tumble" , based on the rotation 
of its flagellae. Technically, it would be possible to run the detailed model for 
all space and time, and observe the macroscopic variables of interest, but this 
would be prohibitively expensive. It is known, however, that, under certain 
conditions, one can write a closed deterministic model for the evolution of the 
concentration of the bacteria as a function of space and time. 

The recently proposed gap-tooth scheme |KGK02] can then be used instead 
of performing stochastic time integration in the whole domain. A number of 
small intervals (boxes, teeth), separated by large gaps, are introduced; they 
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qualitatively correspond to mesh points for a traditional, continuum solution 
of the unavailable chemotaxis equation. The scheme works as follows: We 
first choose a number of macroscopic grid points. We then choose a small 
interval around each grid point; initialize the fine scale, microsopic solver within 
each interval consistently with the macroscopic initial conditions; and provide 
each box with appropriate (as we will see, to some extent artificial) boundary 
conditions. Subsequently, we use the microscopic model in each interval to 
simulate evolution until time At, and obtain macroscopic information (e.g. by 
computing the average density in each box) at time At. This amounts to a 
coarse time-At map; this procedure is then repeated. 

The generalized Godunov scheme of E and Engquist |EE03| also solves an 
unavailable macroscopic equation by repeated calls to a microscopic code; how- 
ever, the assumption is made that the unavailable equation can be written in 
conservation form. In the gap-tooth scheme discussed here, the microscopic 
computations are performed without assuming such a form for the "right-hand- 
side" of the unavailable macroscopic equation; we evolve the detailed model 
in a subset of the domain, and try to recover macroscopic evolution through 
interpolation in space and extrapolation in time. 

We have showed analytically, in the context of numerical homogenization, 
that the gap-tooth scheme is close to a finite difference scheme for the homoge- 
nized equation [2^R03 . However, that analysis employed simulations using an 
algebraic constraint, ensuring that the initial macroscopic gradient is preserved 
at the boundary of each box over the time-step At. This requires altering an 
existing microscopic code, so as to impose this macroscopically-inspired con- 
straint. This may be impractical (e.g. if the macroscopic gradient has to be 
estimated), undesirable (e.g. if the development of the code is expensive and 
time-consuming) or even impossible (e.g. if the microscopic code is a legacy 
code). Generally, a given microscopic code allows us to run with a set of pre- 
defined boundary conditions. It is highly non-trivial to impose macroscopically 
inspired boundary conditions on such microscopic codes, see e.g. |LLY98| for a 
control-based strategy. This can be circumvented by introducing buffer regions 
at the boundary of each small box, which shield the short-time dynamics within 
the computational domain of interest from boundary effects. One then uses the 
microscopic code with its built-in boundary conditions. 

Here, we show we can study the gap-tooth scheme (with buffers) through 
its numerically obtained damping factors, by estimating its eigenvalues. Inte- 
gration with nearby coarse initial conditions is used to estimate matrix-vector 
products of the linearization of the coarse time- At map with known perturbation 
vectors; these are integrated in matrix-free iterative methods such as Arnoldi 
eigensolvers. For a standard diffusion problem, we show that the eigenvalues of 
the gap-tooth scheme are approximately the same as those of the finite differ- 
ence scheme. When we impose Dirichlet boundary conditions at the boundary 
of the buffers, we show that the scheme converges to the standard gap-tooth 
scheme for increasing buffer size. 
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2 Physical Problem/ Governing Equations 



We consider a general reaction-convection-difFusion equation with a dependence 
on a small parameter e, 

— i) = / ( u(x, t), -q^u{x, t), q^H^i t),x,-j , (1) 

with initial condition u{x, 0) = ^0(2;) and Dirichlet boundary conditions u{0, t) = 
Vo and u{l,t) = vi. We further assume that / is 1-periodic in y = |. 

Since we are only interested in the macroscopic (averaged) behavior, let us 
define an averaging operator for u{x,t) as follows 

U{x,t) := Sh{u){x,t) = %i(e,t)de (2) 

This operator replaces the unknown function with its local average in a small 
box of size h » e around each point. If h is sufficiently small, this amounts to 
the removal of the microscopic oscillations of the solution, retaining its macro- 
scopically varying components. 

The averaged solution U{x,t) satisfies an (unknown) macroscopic partial 
differential equation, 

§-U{x, t)=F (u{x, t), ^U{x, t), ^U{x, t), X- h)j , (3) 

which does not depend on the small scale, but instead has a dependence on the 
box width h. 

The goal of the gap-tooth scheme is to approximate the solution U{x,t), 
while only making use of the detailed model Q). For illustration purposes, 
consider as a microscopic model the constant coefficient diffusion equation, 

d 92 

—u{x,t) = a*—u{x,t), (4) 

In this case both U{x,t) and u{x,t) satisfy l@J. The microscopic and macro- 
scopic models are the same, which allows us to focus completely on the method 
and its implementation. 



3 Multiscale/Multiresolution Method 
3.1 The gap-tooth scheme 

Suppose we want to obtain the solution of the unknown equation Q on the 
interval [0,1], using an equidistant, macroscopic mesh n(Aa;) := {0 = xq < 
xi = xq + Ax < . . . < xn = 1}. To this end, consider a small interval {tooth, 
box) of length h « Ax centered around each mesh point, and let us perform a 
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time integration using the microscopic model (QJ in each box. We provide each 
box with boundary conditions and initial condition as follows. 

Boundary conditions. Since the microscopic model ^ is diffusive, it 
makes sense to impose a fixed gradient at the boundary of each small box for 
a time At for the macroscopic function U{x,t). The value of this gradient is 
determined by an approximation of the concentration profile by a polynomial, 
based on the (given) box averages U^, i = 1, . . . ,N. 

^(x, ^7^) ~ Pj^ (x, ^n): X G [Xi — , X^ -|- — ], 

where p^{x\tn) denotes a polynomial of (even) degree k. We require that the 
approximating polynomial has the same box averages in box i and in | boxes 
to the left and to the right. This gives us 



1 k k 

-hi ^^(^=^")de = f^. ■^--2'--'2 



One can easily check that 

k_ k 
5,(pf)(x,i„) = UIX.Ll^ix), Ll^ix) = n (6) 

^\ [Xi+J -x,+i) 

where Lfj (x) denotes a Lagrange polynomial of degree k. The derivative of 
Pi{x, tn) at Xi ± ^ is subsequently used as a Neumann boundary condition. 



^P?(x;t„) 



(7) 



In |SKR03j , it is shown how to enforce the macroscopic gradient to be constant 
when the microscopic model exhibits fast oscillations. 

Initial condition. For the time integration, we must impose an initial 
condition u^{x, i„) in each box [x, — -I, Xi + -1], at time We require ■5,*(x, i„) 
to satisfy the boundary condition and the given box average. We choose a 
quadratic polynomial, centered around the coarse mesh point, 

u\x, In) = a(x — Xi)^ + b{x - Xi) + c. (8) 
Using the constraints {Tj) and requiring j- J'^'^h 'u^{(,itn)d^ = U^, we obtain 

a 2h ' 2 ' ^ 24^ ' ' ^ ' 

The algorithm. The complete gap-tooth algorithm to proceed from tn to 
tn+i = tn + At is given below: 
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1. At time construct the initial condition u'(x,t„), i = 0, . . . , iV, using 
the box averages U" {j — 0, . . . , N) as defined in 

2. Compute u^{x, t) by solving the equation JQ) imtil time tn+i = t + At with 
Neumann boundary conditions 0. 

3. Compute the box average ?7"^^ at time tn+i- 

It is clear that this amounts to a "coarse to coarse" time-At map. We write 
this map as follows, 

=5fe(;7";t„ + At), (10) 

where S represents the numerical time-stepping scheme for the macroscopic 
(coarse) variables and k denotes the degree of interpolation. 

We emphasize that the scheme is also applicable if the microscopic model 
is not a partial differential equation. In this case, we replace step 2 with a 
coarse time-stepper, based on the lift-run-restrict procedure that was outlined in 
[CKT02 . Numerical experiments using this algorithm are presented in GLK03. 



IGK02| . Figure 1(a) gives a schematic representation of the algorithm. 



3.2 The gap-tooth scheme with buffers 

We already mentioned that, in many cases, it is not possible or convenient to 
constrain the macroscopic gradient. However, the only crucial issue is that the 
detailed system in each box should evolve as if it were embedded in a larger do- 
main. This can be effectively accomplished by introducing a larger box of size 
H >> h around each macroscopic mesh point, but still only use (for macro- 
purposes) the evolution over the smaller, "inner" box. This is illustrated in 



figure 1(b) Lifting and evolution (using arbitrary outer boundary conditions) 
are performed in the larger box; yet the restriction is done by taking the average 
of the solution over the inner, small box. The goal of the additional computa- 
tional domains, the buffers, is to buffer the solution inside the small box from 
outer boundary effects. This can be accomplished over short enough times, 
provided the buffers are large enough; analyzing the method is tantamount to 
making these statements quantitative. 

The idea of using a buffer region was also used in the multi-scale finite 
element method (oversampling) of Hon jHW97| to eliminate the boundary layer 
effect; also Hadjiconstantinou makes use of overlap regions to couple a particle 
simulator with a continuum code |Had99| . If the microscopic code allows a 
choice of different types of "outer" microscopic boundary conditions, selecting 
the size of the buffer may also depend on this choice. 



4 Results 

We first show analytically and numerically that the standard gap-tooth scheme 
converges for equation Q). We then analyze convergence of the scheme with 
buffers and Dirichlct boundary conditions through its damping factors. 
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(a) The gap-tooth scheme 



(b) Introduction of buffers 



Figure 1: A schematic representation of the two variants of the gap-tooth scheme 
(with and without buffer boxes). 

4.1 Convergence of the gap-tooth scheme 

Theorem 1. The gap-tooth scheme, applied to equation ^ with exact, an- 
alytical integration within the boxes, and boundary conditions defined through 
interpolating polynomials of (even) order k, is equivalent to a finite difference 
discretization of equation ^ of order k central differences in space and an ex- 
plicit Euler time step. 

Proof. When using exact (analytic) integration in each box, we can find an 
exphcit formula for the gap-tooth time-stepper. The initial profile is given by 
Ui(x,tn) — a{x — + b{x — Xi) + c. Due to the Neumann bomidary 
conditions, time integration can be done analytically, using 

u{x, tn + At) = a{x — Xif' -f- bix — Xi) + c + 2a • a* At 

Averaging this profile over the box gives the following time-stepper for the box 
averages, 

t/« = U^ + a* ^' At. 

h 

We know that sf — j^Pi{x,tn)\,_^,_j_h , where p'^{x,tn) is determined by ©. 
One can easily verify that 

^5.(,?)(.,t„).£t^, 

is a fc-th order approximation of which concludes the proof. □ □ 

As an example, we apply the gap-tooth scheme to the diffusion equation 
(0} with a* ~ 1. We choose an initial condition U{x,t) 1 — \2x — 1\, with 
Dirichlet boundary conditions, and show the result of a fourth-order gap-tooth 
simulation with Ax = 0.05, At = 5 • 10"'^ and h ~ 0.01. Inside each box, we 
used a second order finite difference scheme with microscopic spatial mesh size 
6x^1- 10^3 and 5t ^ b ■ IQ-'' . The resuhs are shown in figure El 
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(a) Numerical solution (b) Difference with FD 

Figure 2: The gap-tooth scheme of fourth order for eq. at i = 0, 4-10^'^, . . . , 2- 
10-2. 

4.2 Damping factors 

Convergence results are typically established by proving consistency and stabil- 
ity. If one can prove that the error in each time step can be made arbitrarily 
small by refining the spatial and temporal mesh size, and that an error made at 
time tn does not get amplified in future time-steps, one has proved convergence. 
This requires the solution operator to be stable as well. 

In the abscense of explicit formulas, one can examine the damping factors 
of the time-stepper. If, for decreasing mesh sizes, all (finitely many) eigenvalues 
and eigenfunctions of the time-stepper converge to the dominant eigenvalues 
and eigenfunctions of the time evolution operator, one expects the solution of 
the scheme to converge to the true solution of the evolution problem. 

Consider equation with Dirichlet boundary conditions u{0,t) — and 
u{l,t), and denote its solution at time t by the time evolution operator 

u{x,t) = s{uo{x);t), (11) 

We know that 

s(sin(TO7ra;); t) = e-^™'^)^* sin(m7ra;), m £ N. 

Therefore, if we consider the time evolution operator over a fixed time t, s{-, t), 
then this operator has eigenftmctions sin(m7rx), with resp. eigenvalues 

A™ = e-('"'^)'*; (12) 

A good (finite difference) scheme approximates well all eigenvalues whose eigen- 
functions can be represented on the given mesh. We note that it is possible to 
decouple the time horizon t from the gap-tooth (or finite-difference) time-step 
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Figure 3: Comparison between the damping factors (left) and the eigenfunction 
corresponding to eigenvalue A3 (right) of the exact solution (full line), the finite 
difference approximation (dashed) and the gap-tooth scheme (dotted). The 
eigenfunction of the gap-tooth scheme is indistinguishable of the finite difference 
eigenfunction. 

At, in order to study the effect of different discretizations on the same reporting 
horizon. 

Since the operator defined in Hll|l is linear, the numerical time integration is 
equivalent to a matrix-vector product. Therefore, we can compute the eigenval- 
ues using matrix-free linear algebra techniques, even for the gap-tooth scheme, 
for which it might not even be possible to obtain a closed expression for the 
matrix. We note that this analysis gives us an indication about the quality of 
the scheme, but it is by no means a proof of convergence. 

We illustrate this with the computation of the eigenvalues of the gap-tooth 
scheme with Neumann box boundary conditions. In this case, we know from 
theorem 1 that these eigenvalues should correspond to the eigenvalues of a finite 
difference scheme on the same mesh. We compare the eigenvalues of the gap- 
tooth scheme of order k — 2 for equation with diffusion coefficient a* = 
0.45825686. As method parameters, we choose Ax — 0.05, h = 5 ■ 10^^, At = 
2.5- 10~* for a time horizon t = 4- 10^"^, which corresponds to 16 gap-tooth steps. 
Inside each box, we use a finite difference scheme of order 2 with 6x = 1-10~^ and 
an implicit Euler time-step of 5 ■ 10"^. We compare these eigenvalues to those 
the finite difference scheme with Ax = 0.05 and At = 2.5 ■ 10"'', and with the 
dominant eigenvalues of the "exact" solution (a finite difference approximation 
with Ax = 1- 10^3 and At = 1 • 10"'^). The resuh is shown in figure El The 
difference between the finite difference approximation and the gap-tooth scheme 
in the higher modes, which should be zero according to theorem 1, is due to the 
numerical solution inside each box and the use of numerical quadrature for the 
average. 

We now examine the effect of introducing a buffer region, as described in 
section IX^ We consider again equation Q with a* — 0.45825686, and we take 
the gap-tooth scheme with parameters Ax ~ 0.05, h ~ 5 ■ 10^^, At — 2.5 • 10~^ 
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Figure 4: Comparison between the damping factors (left) and the eigenfunction 
corresponding to the eigenvalue A3 (right) of the exact solution (full line), the 
finite difference scheme (dashed) and the gap-tooth scheme with buffers (dash- 
dotted lines) for increasing buffer sizes H = 2 ■ 10^^, 3 • fO^^ . . . , f ■ fO^^. 



for a time horizon t = A ■ 10~^, and an internal time-stepper as above. We 
introduce a buffer region of size H, and we impose Dirichlet boundary conditions 
at the outer boundary of the buffer region. Lifting is done in identically the 
same way as for the gap-tooth scheme without buffers; we only use Q as the 
initial condition in the larger box [x^ — -^i + ■§■]■ We compare the eigenvalues 
again with the equivalent finite difference scheme and the exact solution, for 
increasing sizes of the buffer box H. Figure^ shows that, as H increases, the 
eigenvalues of the scheme converge to those of the original gap-tooth scheme. 
We see that, in this case, we would need a buffer of size H = 4 ■ fO~^, i.e. 
80% of the original domain, for a good approximation of the damping factors. 
Note that it is conceptually possible, though inefficient, to let the buffer boxes 
overlap. 



5 Summary/Conclusions 

We described the gap-tooth scheme for the numerical simulation of multi-scale 
problems. This scheme simulates the macroscopic behaviour over a macroscopic 
domain when only a microscopic model is explicitly available. In the case of 
diffusion, we showed equivalence of our scheme to standard finite differences of 
arbitrary (even) order, both theoretically and numerically. 

We showed that it is possible, even without analytic formulas, to study the 
properties of the gap-tooth scheme and generalizations through the damping 
factors of the resulting coarse time-At map. We illustrated this for the original 
gap-tooth scheme and for an implementation using Dirichlet boundary condi- 
tions in a buffer box. We showed that, as long as the buffer region is "large 
enough" to shield the internal region from the boundary effects over a time At, 
we get a convergent scheme. Therefore, we are able to use microscopic codes in 
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the gap-tooth scheme without modification. 

In a forthcoming paper, we will explore, using these damping factors for 
many different types of boundary conditions, the relation between the quality 
of the boundary conditions, the size of the buffer region and the time-step be- 
fore reinitialization. We will investigate the trade-off between the effort required 
to impose a particular type of boundary conditions (and the eventual macro- 
scopically inspired control-based strategy) and the efficiency gain due to smaller 
buffer sizes and/or longer possible time-steps before reinitialization. Here, we 
showed that this investigation is made possible by studying the damping factors 
of the resulting coarse time-At map. 
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